#!/usr/local/bin/python

import os
import shutil
import sys
import subprocess
import chipseq_pipeline
from chipseq_pipeline import *
import redmine

# ----------------------------[ Settings ]--------------------------------------

# Setting for analysis
WHITE_LIST= ''  ## Regions to *include* in the filtered BAMs. Set None or '' if no such white list exists and/or every region is to be included.
GQUADS_DB= '/data01/sblab/users/berald01/reference_seqs/gquads_mmusculus_37.bed'
IGV_GENOME= 'mm9'        ## Used to produced .tdf files from bedgraph. See ~/applications/igv/IGVTools/genomes/ for available ids (ids= id w/o .genome. E.g. mm9 = mm9.genome)
HOMER_GENOME= IGV_GENOME  ## homer genome id might be the same as igv (ok: hg18, hg19, mm9, mm8)
ANNOVAR_BUILDVER= IGV_GENOME ## Passed to annotate_variation.pl --buildver
ANNOVAR_DB= '/home/berald01/applications/annovar/annovar/mousedb' ## Database dir for annotate_variation.pl
BEDTOOLS_GENOME= '/home/berald01/applications/bedtools/BEDTools-Version-2.14.3/genomes/mouse.mm9.genome'
FUNCTIONAL_ANNOTATIONS= [GQUADS_DB, '/data01/sblab/users/berald01/reference_seqs/cpgIslandsExt.NCBI37_mm9.bed'] ## List of annotation files (bed) for which to compute read coverage 

# Settings for export results
PGCONN= {'host':'xxx.xx.xxx.xxx', 'database':'sblab', 'user':'me', 'scp_user':'otherme'}
PROJECT_ID= '20111212_chipseq_formylseq_raiberea'   ## This id must be present in table projects
REDMINE_WIKI= 'http://xxx.xx.xxx.xxx:3000/projects/chipseq-formylseq-raiberea-20111118/wiki/Wiki-chipseq-pipeline' ## This page *must* exist and will be __overwritten__.

# These settings probably don't need changing
PGSCHEMA= 'chipseq_pipeline'
DESIGN= 'design.txt'
DESCRIPTION= 'Output from chipseq_pipeline' ## General description for the output files. This string passed to "project_files -d"

## Do not change
DIR_CLEAN= 'bam_clean' ## Where filtered bam files will be (see bam_clean()) 
MACS_PEAK_DIR= 'macs/peaks'
PACKAGE_REP= '/home/berald01/python-packages/' ## Where the chipseq_pipeline package is
# ----------------------------[ Clean-up ]-------------------------------------

if False:
    shutil.rmtree(DIR_CLEAN, ignore_errors= True) ## Output dir from clean_bam()
    shutil.rmtree('homer', ignore_errors= True)
    shutil.rmtree('macs', ignore_errors= True)
    shutil.rmtree('stderr', ignore_errors= True)
    shutil.rmtree('stdout', ignore_errors= True)
    shutil.rmtree('redmine_wiki', ignore_errors= True)
    shutil.rmtree('fastqc', ignore_errors= True)
    shutil.rmtree('bedgraph', ignore_errors= True)
    try: os.remove('commands.log')
    except: pass
    try: os.remove('stderr.log')
    except: pass
    try: os.remove('stdout.log')
    except: pass

# ----------------------------[ Check parameters ]-----------------------------

# ...TODO

# ----------------------------[ Analyses ]-------------------------------------

if False:
    rename_bam('bam_lims', DESIGN)
    fastqc(aln_dir= 'bam_lims', output_dir= 'fastqc/bam_lims', opt= '--threads 8', file_regex= '\.bam$', echo= False)
    clean_bam(bam_lims_dir= 'bam_lims', bam_clean_dir= DIR_CLEAN, file_regex= '\.bam$', output_exts= '.clean.bam', opt= '-q 15 %s' %('-L %s' %WHITE_LIST if WHITE_LIST != '' and WHITE_LIST is not None else '') )
    aln_stats(aln_dir= 'bam_lims', outfile= 'samtools_flagstats.txt', mapq= 15)
    exec_generic(cmd= """ mkdir -p redmine_wiki; ls bam_lims/*.bam | bamqc.py -i - -o bamqc.txt; mv bamqc.txt* redmine_wiki """)
    exec_generic(cmd= """ mkdir -p bedgraph; for bam in `ls bam_clean/*.bam`; do bbam=`basename $(echo $bam)`; genomeCoverageBed -bg -ibam 'bam_clean/'$bbam -g %s > 'bedgraph/'$bbam'.bedgraph'; igvtools tile 'bedgraph/'$bbam'.bedgraph' 'bedgraph/'$bbam'.tdf' %s; rm 'bedgraph/'$bbam'.bedgraph'; done """ %(BEDTOOLS_GENOME, IGV_GENOME), echo= False)
    exec_generic(cmd= """rm bam_clean.readcount.tsv; for b in `ls bam_clean/*.bam`; do bam=`basename $b`; count=`samtools view -c $b`; echo -e "$bam\t$count"; echo -e "$bam\t$count" >> bam_clean.readcount.tsv; done""")
    
if False:
    macs(aln_dir= DIR_CLEAN, design= DESIGN, output_dir= MACS_PEAK_DIR, opt= '--diag', file_regex= '\.bam$', echo= False, output_exts= '.macs', project_id= PROJECT_ID)
    get_macs_negative_peaks(macs_dir= MACS_PEAK_DIR, output_exts= '.macs', output_file= 'macs/macs_negative_peaks.tsv')
    mergePeaks(peak_dir= 'macs/peaks', merged_file= 'macs/mergepeaks/peaks.merged', opt= '-d 500', file_regex= '\.macs_summits\.bed$', echo= False, bed= True)
    exec_generic(cmd= """mkdir -p %(outd)s; cat macs/peaks/*combined.bed | awk -F "\t" '{print $1"\t"$7"\t"$8"\t"$4"\t"$5"\t"$6}' | sort -k 1,1 -k 2,2n -k 3,3n | mergeBed -nms -d -10 -i stdin | awk '{print $1"\t"$2"\t"$3"\tmacs_"$1"_"$2"_"$3"\t"$4"\t+\t%(project_id)s"}' > %(outd)s/macs.mergeBed.bed""" %{'outd': 'macs/mergepeaks', 'project_id': PROJECT_ID}) ## Merge peaks using bedtools/mergeBed. Input comes from combine_macs() using peak starts and ends.
    coverageBed(bam_dir= DIR_CLEAN, ref_bed= 'macs/mergepeaks/macs.mergeBed.bed', coverage_dir= 'macs/coverage', file_regex= '\.bam$', echo= False)
    coverage_matrix(coverage_dir= 'macs/coverage', matrix_file= 'macs/coverage/coverage_matrix.txt', opt= "-x -4 -e '\.bed$' -s '\.coverage\.bed$'", echo= False)
    call_matrix(call_file= 'macs/mergepeaks/macs.mergeBed.bed', matrix_file= None, calls_sep= ';', row_header= 3, calls_index= 4, strip_dir= True, strip_exts= True, sub_regex= '_peak_.*')
    exec_generic(cmd= """mkdir -p macs/annotatepeaks/; annotatePeaks.pl %s %s | tail -n +2 | homerPeaks2bed.py | closestBed -t first -D ref -a stdin -b %s > %s""" %('macs/mergepeaks/macs.mergeBed.bed', HOMER_GENOME, GQUADS_DB, 'macs/annotatepeaks/macs.mergeBed.annotated'))
    paste_mergefiles(mergef= 'macs/annotatepeaks/macs.mergeBed.annotated', covf= 'macs/coverage/coverage_matrix.txt', callf= 'macs/mergepeaks/macs.mergeBed.bed.call_matrix', outfile= 'macs/macs.mergeBed.pasted.bed')    
    exec_generic(cmd= 'mkdir -p macs/gquads; closestBed -d -a %s -b %s > %s' %('macs/annotatepeaks/peaks.annotated.bed', GQUADS_DB, 'macs/gquads/peaks.annotated.quads.bed'), echo= False)

if False:
    " annovar annotation "
    exec_generic(cmd= """awk -F "\\t" '{print $1"\\t"int(($2 + $3)/2)"\\t"int(($2 + $3)/2)"\\t0\\t0\\t"$4}' %(inbed)s > %(inannovar)s; annotate_variation.pl --buildver %(buildver)s  -geneanno  %(inannovar)s %(annovar_db)s """ %{'inbed':'macs/mergepeaks/macs.mergeBed.bed', 'buildver': ANNOVAR_BUILDVER, 'inannovar':'macs/mergepeaks/macs.mergeBed.annovar', 'annovar_db': ANNOVAR_DB})
    exec_generic(cmd= """R CMD BATCH %s/chipseq-pipeline/annovar_variant_function.R """ %(PACKAGE_REP,))

# ---------------------------[ Functional coverage ]---------------------------
if False:
    for bed in FUNCTIONAL_ANNOTATIONS:
        coverageBed(bam_dir= DIR_CLEAN, ref_bed= bed, coverage_dir= os.path.join('functional_coverage', os.path.splitext(os.path.basename(bed))[0]), file_regex= '\.bam$')


if False:
    "  ---------------  findPeaks module is currently deprecated ------------- "
    
    homer_maketagdirs(input_dir= DIR_CLEAN, homertag_dir= 'homer/tag_directories', echo= False, file_regex= '\.bam$')
    homer_findPeaks(tag_dir= 'homer/tag_directories', design= DESIGN, findpeaks_dir= 'homer/findpeaks', opt= '-center', echo= False)
    get_findPeaks_summary('homer/findpeaks', file_regex= '\.findPeaks$', summary_file= 'homer/peak_summary/findpeaks_summary.txt')
    mergePeaks(peak_dir= 'homer/findpeaks', merged_file= 'homer/mergepeaks/peaks.merged', opt= '-d 500', file_regex= '\.findPeaks$', echo= False, bed= False)
    annotatePeaks(peak_file= 'homer/mergepeaks/peaks.merged', annotated_file= 'homer/annotatepeaks/peaks.annotated', opt= '', refgenome= HOMER_GENOME, echo= False)
    coverageBed(bed_dir= 'bed_clean', ref_bed= 'homer/annotatepeaks/peaks.annotated.bed', coverage_dir= 'homer/coverage', file_regex= '\.bed$', echo= False)
    coverage_matrix(coverage_dir= 'homer/coverage', matrix_file= 'homer/coverage/coverage_matrix.txt', opt= "-x -4 -e '\.bed$'", echo= False)
    call_matrix(call_file= 'homer/annotatepeaks/peaks.annotated.bed', matrix_file= None, row_header= 3, calls_index= 6, strip_dir= True, strip_exts= True)
    exec_generic(cmd= 'mkdir -p homer/gquads; closestBed -d -a %s -b %s > %s' %('homer/annotatepeaks/peaks.annotated.bed', GQUADS_DB, 'homer/gquads/peaks.annotated.quads.bed'), echo= False)
        
# ------------------------------[ Postgres Upload ]----------------------------

if False:
    exec_generic("""project_files.py -v -p %s -d "%s" -c "%s" -r -f .""" %(PROJECT_ID, DESCRIPTION,  'host=' +  PGCONN['host'] + ' user=' + PGCONN['user'] + ' dbname=' + PGCONN['database']))

## This part deprectaed
#if False:
#    pg_read_table(infile= DESIGN, schema= PGSCHEMA, table= 'chipseq_pipeline_design', conn= PGCONN, read_table_opt= "sep: '\t', overwrite: True, header: True, comment_char: '#'")
#    pg_read_table(infile= 'homer/peak_summary/findpeaks_summary.txt', schema= PGSCHEMA, table= 'homer_findpeaks_summary', conn= PGCONN, read_table_opt= "sep: '\t', overwrite: True, header: True, comment_char: '#'")
#    pg_upload_annotatePeaks_gquads_bed(infile= 'homer/gquads/peaks.annotated.quads.bed', conn= PGCONN, schema= PGSCHEMA, table= 'homer_peaks_annotated', overwrite= True)
#    pg_upload_annotatePeaks_gquads_bed(infile= 'macs/gquads/peaks.annotated.quads.bed', conn= PGCONN, schema= PGSCHEMA, table= 'macs_peaks_annotated', overwrite= True)
#    pg_read_table(infile= 'redmine_wiki/macs.bed', schema= PGSCHEMA, table= 'macs_peaks', conn= PGCONN, read_table_opt= " overwrite: False, skip: 1, append: True")

# ----------------------------[ Redmine ]--------------------------------------
if True:
    redmine.remove_attachments(redmine_wiki=  REDMINE_WIKI)
    wiki_text="""
_This page and its file(s) attached have been produced by the chipseq pipeline_

%(svn_info)s

{{toc}}

For standard output/error see attached attachment:stdout.log and attachment:stderr.log

h1. Chipseq design

Experimental design as extracted from the design file.

%(design)s

h1. Alignment statistics

Counts computed by @samtools view -c@

%(aln_stats)s

Duplication level (percentage)

%(dup_level)s

See also the attached quality control from FastQC attachment:bam_lims_fastqc.tar.gz

h1. Alignment statistics from bamqc.py

%(bamqc_txt)s

h1. Peak calling: HOMER

For the homer suite see "here":http://biowhat.ucsd.edu/homer/chipseq/

%(homer_IP_summary)s

See also complete summary output attachment:findpeaks_summary.txt.

h1. Peak calling: MACS

For macs see "homepage":http://liulab.dfci.harvard.edu/MACS/ and "documentation":http://liulab.dfci.harvard.edu/MACS/README.html

%(macs_IP_summary)s

See also the complete summary output attachment:macs_xls_long_summary.txt.

%(macs_neg_peaks)s

Number of peaks detected in the _input control_ file using the _treatment_ file as control


!%(macs_hist)s!

!%(macs_boxplot)s!

h2. Annotation plots

These plots classify peaks (peak centers) according to which functional category their region belongs to.
Produced by @annovar@ tool @annotate_variation.pl@. See also "annovar":http://www.openbioinformatics.org/annovar/ for documentation

!%(annovar_pie)s!
!%(annovar_barplot)s!


h1. Log of the executed commands

<pre>
%(cmd_log)s
</pre>
    """ %{'svn_info':         get_svn_info(chipseq_pipeline),
          'design':           redmine.file2table(DESIGN),
          'aln_stats':        redmine.file2table('redmine_wiki/samtools_flagstats.txt'),
          'dup_level':        redmine.file2table('redmine_wiki/fastqc_dupl_lev.txt'),
          'bamqc_txt':        redmine.file2table('redmine_wiki/bamqc.txt', table_title= ''),
          'homer_IP_summary': redmine.file2table('homer/peak_summary/findpeaks_summary.txt.reduced', table_title= 'h3. Summary extracted from findPeaks\n'),
          'macs_IP_summary':  redmine.file2table('redmine_wiki/macs_xls_short_summary.txt', table_title= 'h3. Summary extracted from macs\n'),
          'macs_neg_peaks':   redmine.file2table('macs/macs_negative_peaks.tsv', table_title= 'h3. Number of peaks negative\n'),
          'annovar_pie':      redmine.attach_wiki(redmine_page= REDMINE_WIKI, attachment= 'macs/pie_annotation.png'),
          'annovar_barplot':  redmine.attach_wiki(redmine_page= REDMINE_WIKI, attachment= 'macs/barplot_annotation.png'),
          'macs_hist':        redmine.attach_wiki(redmine_page= REDMINE_WIKI, attachment= 'redmine_wiki/macs_peak_length_trim_hist.png'),
          'macs_boxplot':     redmine.attach_wiki(redmine_page= REDMINE_WIKI, attachment= 'redmine_wiki/macs_peak_summit_pileup_boxplot.png'),
          'tree': 'NA',
          'cmd_log':          redmine.file2text('commands.log')}
    
    redmine.write_wiki(text= wiki_text, redmine_page= REDMINE_WIKI)
    
    redmine.attach_wiki(redmine_page= REDMINE_WIKI, attachment= 'homer/peak_summary/findpeaks_summary.txt')
    redmine.attach_wiki(redmine_page= REDMINE_WIKI, attachment= 'redmine_wiki/samtools_flagstats.txt')
    redmine.attach_wiki(redmine_page= REDMINE_WIKI, attachment= 'redmine_wiki/bam_lims_fastqc.tar.gz')
    redmine.attach_wiki(redmine_page= REDMINE_WIKI, attachment= 'macs/peaks/macs.combined.bed.tar.gz')
    redmine.attach_wiki(redmine_page= REDMINE_WIKI, attachment= 'redmine_wiki/macs_xls_long_summary.txt')
    redmine.attach_wiki(redmine_page= REDMINE_WIKI, attachment= 'stdout.log')
    redmine.attach_wiki(redmine_page= REDMINE_WIKI, attachment= 'stderr.log')

sys.exit()